/* =============================================================================
This program displays the estimated densities of the wage after changes
in the minimum wage for Canadian provinces. This version includes the new
results including a 12 month lag of the minimum wage.

The do-file loads a set of estimates and uses them to generate densities based
on synthetic data (including all categories used to estimate the model) or 
observed data

Written by: James Townsend (latest version November 2014)
Modified by: Oscar Becerra
Modified by: Iain Snoddy
Date: August 2015
============================================================================= */

global dir_path = "C:\Users\jimt2\Dropbox\min wage distr\JOLE revision\RESULTS_in_paper"
global j=1 /*j=1 all, j=2 joiners, j=4 stayers */

cap program drop trial

prog trial, eclass
args filename sheet
ereturn clear
preserve
import excel "${dir_path}/`filename'.xlsx", describe

import excel  "${dir_path}/`filename'.xlsx", clear  firstrow sheet(`sheet')
sort varname
mkmat b, matrix(A)
*matrix list A
levels varname,local(names)
matname A `names',rows(.) explicit


matrix b=A'

ereturn post b , depname(fweight)

ereturn display
restore


end


// General settings

clear
discard
cap log close

clear
set more off
set matsize 800
set scheme s1color
cap set autotabgraphs on


*sysdir set PERSONAL P:\AA_RESEARCHER_TRAINING_INFORMATION\STATA_ADO_PACKAGES
*global workdatadir "P:\Townsend_5609_wages\November2019\WorkingData"
global logdir "${dir_path}\fitted_densities\logfiles"

global plotdata 	"$dir_path\fitted_densities"
*local wkdate = 20170909
local wkdate = 20241213
log using "$logdir\fitted_densities_compare_`wkdate'.txt", replace text


// Define flags:
// 1. "pub" is a flag for testing the code using the public use files and should
//    be always set to zero when using the data in the RDC.
// 2. "switch" = 0 : use the data set to construct fitted densities from the data set.
//		   	   = 1 : use types specified below to fit densities.
local pub    = 0
local switch = 1

// Additional settings
// ==================== //

// Directories for the public and RDC datasets


// Data set for fitted profiles for "types" - alternative is to add data restrictions.
// Definition of a `type'. The following dataset compares different types of workers.
// The types are defined based on the variables included in the duration model. In
// the example below, we compare workers located all along the wage distribution 
// (wagcat), all of them male, 25 years old, less than High School, living in 
// Ontario during November-2009. 

// The do-file compares the same worker given three different scenarios. We have to
// create four types of worker because the timing of the change is different for 
// the leavers specification, so he his own timing type

// Wage categories and types
set obs 820		//	Set observations to 820 (164 wage categories times five types)
gen byte type = irecode(_n,164,328,492,656) + 1
bys type: gen wagcat = _n

// Generate variables determining types (demographics and time)
// ------------------------------------------------------------ //

gen age = 21						// Age

gen educlev = 3						// Education level

gen prov = 35						// Province

gen year = 2008				// Year

gen month = 4						// Month

gen rminw = .						// REAL MINIMUM WAGE!
replace rminw = 6.75 if type == 1	// 1.  
replace rminw = 6.75 if type == 2	// 2. 
replace rminw = 6.75 if type == 3	// 3. 
replace rminw = 7.25 if type == 4   //
replace rminw = 7.25 if type == 5   //


*This is the minimum wage used for the original dynamic terms.
*Note that rminw is (in the backward representation) the current minimum wage
*while rminw is the wage from 2 months ago.
*The cmin* terms produce a forward "scaffold", which vanishes, while the minimum wage
*term min* move to the new minimum wage

gen rminw3 = .						// Current REAL MINIMUM WAGE 
replace rminw3 = 6.75 if type == 1	// 1. 
replace rminw3 = 7.75 if type == 2	// 2. 
replace rminw3 = 8.75 if type == 3	// 3. 
replace rminw3 = 6.75 if type == 4	// 3. 
replace rminw3 = 6.75 if type == 5	// 3. 


gen mwchg = .						// Minimum wage change in past 3 months
replace mwchg = 0 if type == 1
replace mwchg = 0 if type == 2
replace mwchg = 0 if type == 3
replace mwchg = 0 if type == 4
replace mwchg = 0 if type == 5


*These are the terms for the new representation, as well as the 3-12 month effects.
*Note these work differently...the minimum wage effects are always turned on at the 
*current minimum wage, while the minimum wage dynamics are now scaffolds that reach back
*to the original minimum wage (and disappear once the period has elapsed) or reaching
*forward to the new minimum wage. Turning on is achieved by the
*various dummies indicated recent or upcoming changes (i.e. mwchg*)


gen rminw2 = .						// REAL MINIMUM WAGE 1-2 months prior
replace rminw2 = 6.75 if type == 1	// 1. 
replace rminw2 = 6.75 if type == 2	// 2. 
replace rminw2 = 7.25 if type == 3	// 3. 
replace rminw2 = 7.25 if type == 4	// 3. 
replace rminw2 = 7.25 if type == 5	// 3. 


gen mwchg12 = .						// Minimum wage change in past 3 months
replace mwchg12 = 0 if type == 1
replace mwchg12 = 1 if type == 2
replace mwchg12 = 1 if type == 3
replace mwchg12 = 0 if type == 4
replace mwchg12 = 0 if type == 5


gen rminw312 = .						// REAL MINIMUM WAGE 3-12 months prior
replace rminw312 = 6.75 if type == 1	// 1. 
replace rminw312 = 6.75 if type == 2	// 2. 
replace rminw312 = 6.75 if type == 3	// 3. 
replace rminw312 = 6.75 if type == 4	// 3. 
replace rminw312 = 6.75 if type == 5	// 3. 


gen mwchg312 = .						// Minimum wage change in past 3 months
replace mwchg312 = 0 if type == 1
replace mwchg312 = 0 if type == 2
replace mwchg312 = 0 if type == 3
replace mwchg312 = 1 if type == 4
replace mwchg312 = 0 if type == 5


gen rminwf12 = .						// REAL MINIMUM WAGE 1-2 months from now.
replace rminwf12 = 6.75 if type == 1	// 1. 
replace rminwf12 = 7.25 if type == 2	// 2. 
replace rminwf12 = 6.75 if type == 3	// 3. 
replace rminwf12 = 6.75 if type == 4	// 3. 
replace rminwf12 = 6.75 if type == 5	// 3. 


gen mwchgf12 = .						// Minimum wage change in past 3 months
replace mwchgf12 = 0 if type == 1
replace mwchgf12 = 1 if type == 2
replace mwchgf12 = 0 if type == 3
replace mwchgf12 = 0 if type == 4
replace mwchgf12 = 0 if type == 5



gen rminw_lag = .					// REAL MINIMUM WAGE LAGGED 12 MONTHS! - Not used
replace rminw_lag = 10 if type == 1
replace rminw_lag = 10 if type == 2
replace rminw_lag = 11 if type == 3

/***********************************************************From different project
									// Real disabilility benefits
gen disab=.
replace disab = 11992.78 if type == 1
replace disab = 13821.29 if type == 2
replace disab = 12778.59 if type == 3

									// Real SA benefits
gen benefit=.
replace benefit = 5000 if type == 1
replace benefit = 6000 if type == 2
replace benefit = 7000 if type == 3

********************************************************************************/


gen byte dminw = .					// Indicator variable of a change in nominal minimum wage
replace dminw = 0 if type == 1		// For samples 1, 2, and 4 (it means a change in t). For
replace dminw = 0 if type == 2		// leavers, it indicates a change in t + 1
replace dminw = 0 if type == 3


// Based on the variables above, generate variables used in the
// estimation of the duration model
// ------------------------------------------------------------ //

// Generate minimum wage bins. Minimum wages are grouped in 10 cents bins
// from 4.01 to 11.99 dollars. Real minimum wages up to 4 dollars are grouped
// in three categories (up to 3, 3.01 to 3.5 and 3.51 to 4). In practice,
// the real minimum wages are located between XXXXX (bin 18) and XXXX (bin 52)
gen int mincat =   .
replace mincat =   1 if rminw <= 3.0						// Up to 3 dollars
replace mincat =   2 if rminw >  3.0 & rminw <= 3.5			// 3.01 to 3.5 dollars				
replace mincat =   3 if rminw >  3.5 & rminw <= 4.0			// 3.51 to 4.0 dollars				
replace mincat =   4 + (ceil(rminw*10) - 41) if rminw > 4.0	// 4.01(.1)12.0  dollars						
replace mincat =   0 if rminw >= 12							// 12+ dollars

// Test the definition of the bins
table mincat, stat(min rminw) stat(max rminw) stat(count rminw) nformat(%11.6f)

// Generate minimum wage bins. Minimum wages are grouped in 10 cents bins
// from 4.01 to 11.99 dollars. Real minimum wages up to 4 dollars are grouped
// in three categories (up to 3, 3.01 to 3.5 and 3.51 to 4). In practice,
// the real minimum wages are located between XXXX (bin 18) and XXXX (bin 52)
/*************
gen int mincat_lag =   .
replace mincat_lag =   1 if rminw_lag <= 3.0						// Up to 3 dollars
replace mincat_lag =   2 if rminw_lag >  3.0 & rminw_lag <= 3.5		// 3.01 to 3.5 dollars				
replace mincat_lag =   3 if rminw_lag >  3.5 & rminw_lag <= 4.0		// 3.51 to 4.0 dollars				
replace mincat_lag =   4 + (ceil(rminw_lag*10) - 41) if rminw_lag > 4.0	// 4.01(.1)12.0  dollars						
replace mincat_lag =   0 if rminw_lag >= 12							// 12+ dollars

// Test the definition of the bins
table mincat_lag, stat(min rminw_lag) stat(max rminw_lag) stat(count rminw_lag) nformat(%11.6f)
**************/

		gen mincatn=0
		replace mincatn=1 if rminw3<=3
		replace mincatn=2 if rminw3>3 & rminw3<=3.5
		replace mincatn=3 if rminw3>3.5 & rminw3<=4

		local i=1
		while `i'<=76{
			replace mincatn=3+`i' if rminw3>4+(`i'-1)*.1 & rminw3<=4+`i'*.1
			local i=`i'+1
		}

		gen mincat12=0
		replace mincat12=1 if rminw2<=3
		replace mincat12=2 if rminw2>3 & rminw2<=3.5
		replace mincat12=3 if rminw2>3.5 & rminw2<=4

		local i=1
		while `i'<=76{
			replace mincat12=3+`i' if rminw2>4+(`i'-1)*.1 & rminw2<=4+`i'*.1
			local i=`i'+1
		}
		
		
		
		gen mincat312=0
		replace mincat312=1 if rminw312<=3
		replace mincat312=2 if rminw312>3 & rminw312<=3.5
		replace mincat312=3 if rminw312>3.5 & rminw312<=4

		local i=1
		while `i'<=76{
			replace mincat312=3+`i' if rminw312>4+(`i'-1)*.1 & rminw312<=4+`i'*.1
			local i=`i'+1
		}		

gen mincatf12=0
		replace mincatf12=1 if rminwf12<=3
		replace mincatf12=2 if rminwf12>3 & rminwf12<=3.5
		replace mincatf12=3 if rminwf12>3.5 & rminwf12<=4

		local i=1
		while `i'<=76{
			replace mincatf12=3+`i' if rminwf12>4+(`i'-1)*.1 & rminwf12<=4+`i'*.1
			local i=`i'+1
		} 
		
		
sort year month wagcat
		merge m:1 year month wagcat using "$dir_path\dollar_bins"
		drop _merge
    // replace dollar=0
	// replace dollar=0 if type==2
	//	replace dollar=.1 if type==3

// Time. Quarter fixed effects
gen quarter = .
replace quarter = 1 if month >=  1 & month <=  3
replace quarter = 2 if month >=  4 & month <=  6
replace quarter = 3 if month >=  7 & month <=  9
replace quarter = 4 if month >= 10 & month <= 12

label define QUARTER 1 "Jan.-Mar." 2 "Apr.-Jun." 3 "Jul.-Sep." 4 "Oct.-Dec."
label value quarter QUARTER

// Sex

gen byte female = 0

// Age groups
gen byte agecat = .
replace agecat = 1 if age >= 15 & age <= 19
replace agecat = 2 if age >= 20 & age <= 24
replace agecat = 3 if age >= 25 & age <= 34
replace agecat = 4 if age >= 35 & age <= 54
replace agecat = 5 if age >= 55 & age <= 60
drop if agecat == .

label define AGECAT 1 "15-19" 2 "20-24" 3 "25-34" 4 "35-54" 5 "55-60"
label value agecat AGECAT

// Education groups
gen byte edgrp =.
replace edgrp = 1 if educlev >= 0 & educlev <= 2
replace edgrp = 2 if educlev == 3 | educlev == 4
replace edgrp = 3 if educlev >= 5 & educlev <= 7
replace edgrp = 4 if educlev == 8 | educlev == 9
drop if edgrp == .

label define EDGRP 1 "Less than H.S." 2 "H.S. graduate and/or some P.S." 3 "P.S. Diploma or Cert." 4 "B.A. or more"
label value edgrp EDGRP

// Identify sub minimum wage provinces. Generate an indicator variable for 
// Nova Scotia (12) and Ontario (35) for the entire sample, and BC (59) 
// between Nov-2001 and May-2012
gen byte sub = prov == 12 | ///
			   prov == 35 | ///
			  (prov == 59 & inrange(ym(year,month),m(2001m11),m(2012m4)))

// Bin differences (distance) between wage categories and minimum wage 
// categories

/************************************************
JT: Note the following are dated...now wage bins 3-5 and 6-10 bins above minimum wage.
In addition, use 50 cents bins from $1 to $5 above minimum wage.

gen diff     = wagcat - mincat
gen min6b    = diff <  -5				// Wage bin BELOW minwage bin
gen min35b   = diff >= -5 & diff <= -3
gen min12b   = diff == -2 | diff == -1
gen min      = diff ==  0
gen min13a   = diff >=  1 & diff <=  3	// Wage bin ABOVE minwage bin
gen min410a  = diff >=  4 & diff <= 10	// Important - next line has changed definition of min410a
gen min1120a = diff >= 11 & diff <= 20
gen min2130a = diff >= 21 & diff <= 30
gen min3140a = diff >= 31 & diff <= 40
gen min4150a = diff >= 41 & diff <= 50
gen min5160a = diff >= 51 & diff <= 60 

gen byte min6b1    =  dminw == 1 & diff <  -5					// Wage bin BELOW minwage bin
gen byte min35b1   =  dminw == 1 & diff >= -5 & diff <= -3
gen byte min12b1   =  dminw == 1 & diff == -2 | diff == -1
gen byte min1      =  dminw == 1 & diff == 0					// Wage bin AT the minimum wage
gen byte min13a1   =  dminw == 1 & diff >=  1 & diff <=  3		// Wage bin ABOVE minwage bin
gen byte min410a1  =  dminw == 1 & diff >=  4 & diff <= 10		// (Important - this line has changed definition of min410a)
gen byte min1120a1 =  dminw == 1 & diff >= 11 & diff <= 20
gen byte min2130a1 =  dminw == 1 & diff >= 21 & diff <= 30
gen byte min3140a1 =  dminw == 1 & diff >= 31 & diff <= 40
gen byte min4150a1 =  dminw == 1 & diff >= 41 & diff <= 50
gen byte min5160a1 =  dminw == 1 & diff >= 51 & diff <= 60 
******************************************************/


		gen diff=wagcat-mincat /*Difference between wage bin and "current" minimum wage (2 months ago in original specification) */
		gen diff2=wagcat-mincatn /*Difference between wage and new minimum wage */
		gen diff12=wagcat-mincat12 /*Difference between wage and minimum wage 1-2 months ago */
		gen diff312=wagcat-mincat312 
		gen difff12=wagcat-mincatf12 /*Difference between wage and minimum wage 1-2 months ago*/



		gen min6b1=diff<-5
		gen min35b1=diff>=-5&diff<=-3
		gen min12b1=diff==-2|diff==-1
		gen min1=diff==0
		gen min12a1=diff>=1&diff<=2
		gen min35a1=diff>=3&diff<=5
		gen min610a1=diff>=6&diff<=10
		gen min1115a1=diff>=11&diff<=15
		gen min1620a1=diff>=16&diff<=20
		gen min2125a1=diff>=21&diff<=25
		gen min2630a1=diff>=26&diff<=30	
		gen min3135a1=diff>=31&diff<=35
		gen min3640a1=diff>=36&diff<=40	
		gen min4145a1=diff>=41&diff<=45
		gen min4650a1=diff>=46&diff<=50
		*Extended minimum wage effects above $5 follow.
		gen min5155a1=diff>=51&diff<=55
		gen min5660a1=diff>=56&diff<=60	
		gen min6165a1=diff>=61&diff<=65
		gen min6670a1=diff>=66&diff<=70	
		gen min51100a1=diff>=51&diff<=100
		gen min101a1=diff>=101 & diff<=144

		
		
		
		/****IS: Bins around old and new min wage****/
		/****JT: These are the original definitions;
		definitions for 1-2  & 3-12 months since min. wage
		change follow, as well as 1-2 months until new min
		change
		**********************************************/
		gen minold2new=1 if diff>0 & diff2<0
		replace minold2new=0 if minold2new!=1
		gen minnew=diff2==0

		gen min12a2=diff2>=1&diff2<=2
		gen min35a2=diff2>=3&diff2<=5
		gen min610a2=diff2>=6&diff2<=10
		gen min1115a2=diff2>=11&diff2<=15
		gen min1620a2=diff2>=16&diff2<=20
		gen min2125a2=diff2>=21&diff2<=25
		gen min2630a2=diff2>=26&diff2<=30
		gen min3135a2=diff2>=31&diff2<=35
		gen min3640a2=diff2>=36&diff2<=40
		gen min4145a2=diff2>=41&diff2<=45
		gen min4650a2=diff2>=46&diff2<=50
		
		gen cmin6b=mwchg*min6b1
		gen cmin35b=mwchg*min35b1
		gen cmin12b=mwchg*min12b1
		gen cmin=mwchg*min1
		gen cminold2new=mwchg*minold2new
		gen cminnew=mwchg*minnew
		gen cmin12a=mwchg*min12a2
		gen cmin35a=mwchg*min35a2
		gen cmin610a=mwchg*min610a2
		gen cmin1115a=mwchg*min1115a2
		gen cmin1620a=mwchg*min1620a2
		gen cmin2125a=mwchg*min2125a2
		gen cmin2630a=mwchg*min2630a2
		gen cmin3135a=mwchg*min3135a2
		gen cmin3640a=mwchg*min3640a2
		gen cmin4145a=mwchg*min4145a2
		gen cmin4650a=mwchg*min4650a2
		
*Add a "switch here. - new minimum wage variable definitions




		/*gen min6b1=diff<-5
		gen min35b1=diff>=-5&diff<=-3
		gen min12b1=diff==-2|diff==-1
		gen min1=diff==0
		gen min12a1=diff>=1&diff<=2
		gen min35a1=diff>=3&diff<=5
		gen min610a1=diff>=6&diff<=10
		gen min1115a1=diff>=11&diff<=15
		gen min1620a1=diff>=16&diff<=20
		gen min2125a1=diff>=21&diff<=25
		gen min2630a1=diff>=26&diff<=30	
		gen min3135a1=diff>=31&diff<=35
		gen min3640a1=diff>=36&diff<=40	
		gen min4145a1=diff>=41&diff<=45
		gen min4650a1=diff>=46&diff<=50	 */
		
		/****JT: These have been recoded using old vs. current minimum wage ***/
		gen minold2new12=1 if diff<0 & diff12>0
		replace minold2new12=0 if minold2new12!=1
		
		gen minold2new312=1 if diff<0 & diff312>0
		replace minold2new312=0 if minold2new312!=1
		
		*Added - new definition of dynanmics
		gen minold12 = diff12==0
		gen minold312 = diff312==0
		
		gen min6b12=diff12<-5
		gen min35b12=diff12>=-5&diff12<=-3
		gen min12b12=diff12==-2|diff12==-1
		
		gen min6b312=diff312<-5
		gen min35b312=diff312>=-5&diff312<=-3
		gen min12b312=diff312==-312|diff312==-1
		
		
		/****IS: Bins around old and new min wage****/
		gen cmin6b12=mwchg12*min6b12
		gen cmin35b12=mwchg12*min35b12
		gen cmin12b12=mwchg12*min12b12
		gen cmin12=mwchg12*min1  //   JT - now use the current min. wage and work back.
		gen cminold12=mwchg12*minold12
		gen cminold2new12=mwchg12*minold2new12
		gen cminnew12=mwchg12*min1
		gen cmin12a12=mwchg12*min12a1
		gen cmin35a12=mwchg12*min35a1
		gen cmin610a12=mwchg12*min610a1
		gen cmin1115a12=mwchg12*min1115a1
		gen cmin1620a12=mwchg12*min1620a1
		gen cmin2125a12=mwchg12*min2125a1
		gen cmin2630a12=mwchg12*min2630a1
		gen cmin3135a12=mwchg12*min3135a1
		gen cmin3640a12=mwchg12*min3640a1
		gen cmin4145a12=mwchg12*min4145a1
		gen cmin4650a12=mwchg12*min4650a1
		gen cmin5155a12=mwchg12*min5155a1
		gen cmin5660a12=mwchg12*min5660a1
		gen cmin6165a12=mwchg12*min6165a1
		gen cmin6670a12=mwchg12*min6670a1
		
		gen cmin6b312=mwchg312*min6b312
		gen cmin35b312=mwchg312*min35b312
		gen cmin12b312=mwchg312*min12b312
		gen cmin312=mwchg312*min1
		gen cminold312=mwchg312*minold312
		gen cminold2new312=mwchg312*minold2new312
		gen cminnew312=mwchg312*min1	
		gen cmin12a312=mwchg312*min12a1
		gen cmin35a312=mwchg312*min35a1
		gen cmin610a312=mwchg312*min610a1
		gen cmin1115a312=mwchg312*min1115a1
		gen cmin1620a312=mwchg312*min1620a1
		gen cmin2125a312=mwchg312*min2125a1
		gen cmin2630a312=mwchg312*min2630a1
		gen cmin3135a312=mwchg312*min3135a1
		gen cmin3640a312=mwchg312*min3640a1
		gen cmin4145a312=mwchg312*min4145a1
		gen cmin4650a312=mwchg312*min4650a1
		gen cmin5155a312=mwchg312*min5155a1
		gen cmin5660a312=mwchg312*min5660a1
		gen cmin6165a312=mwchg312*min6165a1
		gen cmin6670a312=mwchg312*min6670a1
		
/* Definitions for models with forward looking terms - JT */

		gen minold2newf12= diff>0 & difff12<0

		gen minnewf12=difff12==0

		gen min12a2f=difff12>=1&difff12<=2
		gen min35a2f=difff12>=3&difff12<=5
		gen min610a2f=difff12>=6&difff12<=10
		gen min1115a2f=difff12>=11&difff12<=15
		gen min1620a2f=difff12>=16&difff12<=20
		gen min2125a2f=difff12>=21&difff12<=25
		gen min2630a2f=difff12>=26&difff12<=30
		gen min3135a2f=difff12>=31&difff12<=35
		gen min3640a2f=difff12>=36&difff12<=40
		gen min4145a2f=difff12>=41&difff12<=45
		gen min4650a2f=difff12>=46&difff12<=50
		gen min5155a2f=difff12>=51&difff12<=55
		gen min5660a2f=difff12>=56&difff12<=60
		gen min6165a2f=difff12>=61&difff12<=65
		gen min6670a2f=difff12>=66&difff12<=70
		
	
		gen cmin6bf12=mwchgf12*min6b1
		gen cmin35bf12=mwchgf12*min35b1
		gen cmin12bf12=mwchgf12*min12b1
		gen cminf12=mwchgf12*min1
		gen cminnewf12=mwchgf12*minnewf12
		*gen cminoldf12=mwchgf12*minoldf12
		gen cminold2newf12=mwchgf12*minold2newf12
		gen cmin12af12=mwchgf12*min12a2f
		gen cmin35af12=mwchgf12*min35a2f
		gen cmin610af12=mwchgf12*min610a2f
		gen cmin1115af12=mwchgf12*min1115a2f
		gen cmin1620af12=mwchgf12*min1620a2f
		gen cmin2125af12=mwchgf12*min2125a2f
		gen cmin2630af12=mwchgf12*min2630a2f
		gen cmin3135af12=mwchgf12*min3135a2f
		gen cmin3640af12=mwchgf12*min3640a2f
		gen cmin4145af12=mwchgf12*min4145a2f
		gen cmin4650af12=mwchgf12*min4650a2f
		gen cmin5155af12=mwchgf12*min5155a2f
		gen cmin5660af12=mwchgf12*min5660a2f
		gen cmin6165af12=mwchgf12*min6165a2f
		gen cmin6670af12=mwchgf12*min6670a2f
					

// Split the wage bins in larger groups
gen cut = wagcat


/*************************
gen byte part0 = .
	replace part0 =  1 if cut <=  41
	replace part0 =  2 if cut >=  42 & cut <=  51
	replace part0 =  3 if cut >=  52 & cut <=  61
	replace part0 =  4 if cut >=  62 & cut <=  73
	replace part0 =  5 if cut >=  74 & cut <=  85
	replace part0 =  6 if cut >=  86 & cut <=  98
	replace part0 =  7 if cut >=  99 & cut <= 111
	replace part0 =  8 if cut >=  112 & cut <= 126
	replace part0 =  9 if cut >= 127 & cut <= 142
	replace part0 = 10 if cut >= 143 & cut <= 163
	replace part0 = 11 if cut == 164 
************************/

		/****TL: FIVE SEGMENTS INSTEAD OF TEN****/
			gen part0=.
			  replace part0=1 if cut<=63  /*up to $10*/
			  replace part0=2 if cut>=64&cut<=88 /*$10 to $12.5*/
			  replace part0=3 if cut>=89&cut<=113 /*$12.5 to $15.0*/
			  replace part0=4 if cut>=114&cut<=138 /*$15 to $17.5*/
			  replace part0=5 if cut>=139&cut<=163 /*$17.5 to $20*/
	
	
// Generate variables used in the estimation
local i = 1
while (`i' <= 10) {
	// Province X wage bins (Base category Ontario)
	gen byte p10_`i'= prov == 10 & part0 ==`i' 
	gen byte p11_`i'= prov == 11 & part0 ==`i' 
	gen byte p12_`i'= prov == 12 & part0 ==`i' 
	gen byte p13_`i'= prov == 13 & part0 ==`i' 
	gen byte p24_`i'= prov == 24 & part0 ==`i' 
	gen byte p46_`i'= prov == 46 & part0 ==`i' 
	gen byte p47_`i'= prov == 47 & part0 ==`i' 
	gen byte p48_`i'= prov == 48 & part0 ==`i' 
	gen byte p59_`i'= prov == 59 & part0 ==`i' 

	gen f_`i' = female & part0 ==`i'
	
	// Age groups X wage bins (Base category agegrp == 4)
	gen byte a1_`i'= agecat == 1 & part0 ==`i'
	gen byte a2_`i'= agecat == 2 & part0 ==`i'
	gen byte a3_`i'= agecat == 3 & part0 ==`i'
	gen byte a5_`i'= agecat == 5 & part0 ==`i'

	// Education groups X wage bins (Base category edgrp == 3)
	gen byte e1_`i'= edgrp == 1 & part0==`i'
	gen byte e2_`i'= edgrp == 2 & part0==`i'
	gen byte e4_`i'= edgrp == 4 & part0==`i'
	
	// Quarters X wage bins
	gen byte q1_`i'=quarter==1& part0==`i'
	gen byte q2_`i'=quarter==2& part0==`i'
	gen byte q4_`i'=quarter==4& part0==`i'	

	// Month
/*	gen byte m1_`i' = month == 1 & part0 == `i'
	gen byte m2_`i' = month == 2 & part0 == `i'
	gen byte m3_`i' = month == 3 & part0 == `i'
	gen byte m4_`i' = month == 4 & part0 == `i'
	gen byte m5_`i' = month == 5 & part0 == `i'
	gen byte m6_`i' = month == 6 & part0 == `i'
	gen byte m7_`i' = month == 7 & part0 == `i'
	gen byte m8_`i' = month == 8 & part0 == `i'
	gen byte m9_`i' = month == 9 & part0 == `i'
	gen byte m10_`i' = month == 10 & part0 == `i'
	gen byte m11_`i' = month == 11 & part0 == `i'
	gen byte m12_`i' = month == 12 & part0 == `i'
*/

	// Year X wage bins (Base category year == 1999)
	gen byte y96_`i'= year == 1996 & part0 ==`i' 
	gen byte y97_`i'= year == 1997 & part0 ==`i' 
	gen byte y98_`i'= year == 1998 & part0 ==`i' 
	gen byte y00_`i'= year == 2000 & part0 ==`i'
	gen byte y01_`i'= year == 2001 & part0 ==`i' 
	gen byte y02_`i'= year == 2002 & part0 ==`i' 
	gen byte y03_`i'= year == 2003 & part0 ==`i' 
	gen byte y04_`i'= year == 2004 & part0 ==`i'
	gen byte y05_`i'= year == 2005 & part0 ==`i' 
	gen byte y06_`i'= year == 2006 & part0 ==`i' 
	gen byte y07_`i'= year == 2007 & part0 ==`i' 
	gen byte y08_`i'= year == 2008 & part0 ==`i'
	gen byte y09_`i'= year == 2009 & part0 ==`i' 
	gen byte y10_`i'= year == 2010 & part0 ==`i' 
	gen byte y11_`i'= year == 2011 & part0 ==`i' 
	gen byte y12_`i'= year == 2012 & part0 ==`i'
	gen byte y13_`i'=year==2013 & part0 ==`i'
	gen byte y14_`i'=year==2014 & part0 ==`i'
	gen byte y15_`i'=year==2015 & part0 ==`i'
	gen byte y16_`i'=year==2016 & part0 ==`i'

	/*************************************
	gen ben_`i' = benefit * (part0 ==`i')
	gen disab_`i' = disab * (part0 ==`i')
	*************************************/
	
/*****TL: ADD PROVINCIAL TRENDS****/
		*gen byte ont_`i'=prov==35 & part0 ==`i' 
		gen tr10_`i'=(year-1999)*p10_`i' 
		gen tr11_`i'=(year-1999)*p11_`i' 
		gen tr12_`i'=(year-1999)*p12_`i' 
		gen tr13_`i'=(year-1999)*p13_`i' 
		gen tr24_`i'=(year-1999)*p24_`i' 
		*gen tr35_`i'=(year-1999)*ont_`i' 
		gen tr46_`i'=(year-1999)*p46_`i' 
		gen tr47_`i'=(year-1999)*p47_`i' 
		gen tr48_`i'=(year-1999)*p48_`i' 
		gen tr59_`i'=(year-1999)*p59_`i' 

		gen tq10_`i'=(year-1999)^2*p10_`i' 
		gen tq11_`i'=(year-1999)^2*p11_`i' 
		gen tq12_`i'=(year-1999)^2*p12_`i' 
		gen tq13_`i'=(year-1999)^2*p13_`i' 
		gen tq24_`i'=(year-1999)^2*p24_`i' 
		*gen tq35_`i'=(year-1999)^2*ont_`i' 
		gen tq46_`i'=(year-1999)^2*p46_`i' 
		gen tq47_`i'=(year-1999)^2*p47_`i' 
		gen tq48_`i'=(year-1999)^2*p48_`i' 
		gen tq59_`i'=(year-1999)^2*p59_`i' 

		/******ADD PROV*BIN AND YEAR*BIN CONTROLS*****/
		gen bp10_`i'=(wagcat-(14+25*`i'))*p10_`i' 
		gen bp11_`i'=(wagcat-(14+25*`i'))*p11_`i' 
		gen bp12_`i'=(wagcat-(14+25*`i'))*p12_`i' 
		gen bp13_`i'=(wagcat-(14+25*`i'))*p13_`i' 
		gen bp24_`i'=(wagcat-(14+25*`i'))*p24_`i' 
		gen bp46_`i'=(wagcat-(14+25*`i'))*p46_`i' 
		gen bp47_`i'=(wagcat-(14+25*`i'))*p47_`i' 
		gen bp48_`i'=(wagcat-(14+25*`i'))*p48_`i' 
		gen bp59_`i'=(wagcat-(14+25*`i'))*p59_`i' 

		*gen by96_`i'=(wagcat-(14+25*`i'))*y96_`i'
		gen by97_`i'=(wagcat-(14+25*`i'))*y97_`i' 
		gen by98_`i'=(wagcat-(14+25*`i'))*y98_`i' 
		gen by00_`i'=(wagcat-(14+25*`i'))*y00_`i' 
		gen by01_`i'=(wagcat-(14+25*`i'))*y01_`i' 
		gen by02_`i'=(wagcat-(14+25*`i'))*y02_`i' 
		gen by03_`i'=(wagcat-(14+25*`i'))*y03_`i' 
		gen by04_`i'=(wagcat-(14+25*`i'))*y04_`i' 
		gen by05_`i'=(wagcat-(14+25*`i'))*y05_`i' 
		gen by06_`i'=(wagcat-(14+25*`i'))*y06_`i' 
		gen by07_`i'=(wagcat-(14+25*`i'))*y07_`i' 
		gen by08_`i'=(wagcat-(14+25*`i'))*y08_`i' 
		gen by09_`i'=(wagcat-(14+25*`i'))*y09_`i' 
		gen by10_`i'=(wagcat-(14+25*`i'))*y10_`i' 
		gen by11_`i'=(wagcat-(14+25*`i'))*y11_`i' 
		gen by12_`i'=(wagcat-(14+25*`i'))*y12_`i' 
		gen by13_`i'=(wagcat-(14+25*`i'))*y13_`i' 
		gen by14_`i'=(wagcat-(14+25*`i'))*y14_`i' 
		gen by15_`i'=(wagcat-(14+25*`i'))*y15_`i' 
		gen by16_`i'=(wagcat-(14+25*`i'))*y16_`i' 
	
	
	local i=`i'+1
}

/*****TL: ADD A FULL SET OF prov*year DUMMIES (SEGMENT 1)***/
/*****JT: 20 values per province, starting with NFLD. Ontario is fifth province, so py101 is Ontario in 1997, etc.) ***/
		
		
		
/*************
		gen provyear=100*prov+(year-1997) if part0==1
		replace provyear=9999 if part0>=2
		quietly tab provyear, gen(py)

		/*NOTE: ONTARIO IN 1999 IS py103*/
		/*Note: without 1996, py200 is the last valid dummy (py201 is for obs in the second and higher partition */
*******/

	


		forvalues i = 1/201 {
		gen py`i'=0	
		}
		
		gen p=.
		replace p=0 if prov==10
		replace p=1 if prov==11
		replace p=2 if prov==12
		replace p=3 if prov==13
		replace p=4 if prov==24
		replace p=5 if prov==35
		replace p=6 if prov==46
		replace p=7 if prov==47
		replace p=8 if prov==48
		replace p=9 if prov==59
		
		forvalues p=0/9 {
		forvalues y =1/20 {
		local py=(p*20)+`y'	
		replace py`py' = 1 if p==`p' & year==(1996+`y') 
		
		}
		}
		
	
		*replace py105=1 /*Ontario in 2002*/

		
/*		gen cmin6b=mwchg*min6b1
		gen cmin35b=mwchg*min35b1
		gen cmin12b=mwchg*min12b1
		gen cmin=mwchg*min1
		gen cminold2new=mwchg*minold2new
		gen cminnew=mwchg*minnew
		gen cmin12a=mwchg*min12a2
		gen cmin35a=mwchg*min35a2
		gen cmin610a=mwchg*min610a2
		gen cmin1115a=mwchg*min1115a2
		gen cmin1620a=mwchg*min1620a2
		gen cmin2125a=mwchg*min2125a2
		gen cmin2630a=mwchg*min2630a2
		gen cmin3135a=mwchg*min3135a2
		gen cmin3640a=mwchg*min3640a2
		gen cmin4145a=mwchg*min4145a2
		gen cmin4650a=mwchg*min4650a2 */
		
		gen dolm=dollar*min1
		
/* minimum wage effects interacted with magnitude for models that need them */


/*****JT: Add minimum wage data for interactions*****/
/*		merge m:1 prov year month using $workdatadir\rminw_rminw3
		drop if _merge==2
		drop _merge */
	
*replace rminw3=8.75 if type==2	// Used for proportionality representation, 
*which follows:
	
gen		min1_mw=min1*rminw3
gen  	min6b1_mw=min6b1*rminw3
gen  	min35b1_mw=min35b1*rminw3
gen  	min12b1_mw=min12b1*rminw3
gen  	min12a1_mw=min12a1*rminw3
gen  	min35a1_mw=min35a1*rminw3
gen  	min610a1_mw=min610a1*rminw3
gen  	min1115a1_mw=min1115a1*rminw3  
gen 	min1620a1_mw=min1620a1*rminw3
gen 	min2125a1_mw=min2125a1*rminw3
gen 	min2630a1_mw=min2630a1*rminw3
gen		min3135a1_mw=min3135a1*rminw3
gen 	min3640a1_mw=min3640a1*rminw3
gen 	min4145a1_mw=min4145a1*rminw3
gen  	min4650a1_mw=min4650a1*rminw3
gen	cmin_mw=cmin*rminw3
*could interact everything above here with the new minimum wage (i.e. rminw3)
gen		cminold2new_mw=cminold2new*rminw3
*gen cminold2new_mw3=cminold2new_rminw3 
gen 	cminnew_mw=cminnew*rminw3
gen 	cmin6b_mw=cmin6b*rminw3 
gen		cmin35b_mw=cmin35b*rminw3
gen 	cmin12b_mw=cmin12b*rminw3
gen 	cmin12a_mw=cmin12a*rminw3
gen 	cmin35a_mw=cmin35a*rminw3
gen 	cmin610a_mw=cmin610a*rminw3
gen		cmin1115a_mw=cmin1115a*rminw3 
gen 	cmin1620a_mw=cmin1620a*rminw3
gen		cmin2125a_mw=cmin2125a*rminw3
gen 	cmin2630a_mw=cmin2630a*rminw3
gen 	cmin3135a_mw=cmin3135a*rminw3
gen 	cmin3640a_mw=cmin3640a*rminw3
gen  	cmin4145a_mw=cmin4145a*rminw3 
gen 	cmin4650a_mw=cmin4650a*rminw3		

*New representation minimum wage terms:




// Dummy variables by wage bin
local i=1
while `i'<=164 {
	gen byte bin`i' = cut == `i'
	local i=`i'+1
}

// Remain?
gen remain = 1

sum, sep(0)




save "${plotdata}\synthetic_SA.dta",replace



// -------------------------------------------------------------------------- //
// AFTER THE DATA PROCESSING STAGE, LOAD ESTIMATION RESULTS TO COMPUTE THE 
// DENSITIES FOR THE SELECTED SAMPLE. FOR THIS PART OF THE CODE, SELECT THE
// SPECIFICATION BASED ON SEX AND EARNINGS TYPE. THE FILE LOOPS OVER SUBSAMPLES
// -------------------------------------------------------------------------- //
clear
set more off
global resultdir 	"P:\Townsend_5609_wages\March_2024\results"
global out_plots	"P:\Townsend_5609_wages\March_2024\fitted_densities"
// Names of the files containing the parameter estimates
// The structure of the name of the file is 20150224_`sex'_`k'_`j'_dminw
global ests ${resultdir}/`file'.ster

use "${plotdata}\synthetic_SA.dta",clear

// Define local variables determining the sample/estimation to be used
// local sex = 1	// 1 if male, 2 if female
// local k   = 1	// 1 = full data set; 2 = only paid by the hour; 3 = excluding those earning tips
// local j   = 1	// 1 = full data set; 2 = New hires (tenure < 1yr); 3 = Leavers in period 2;  4 = continuers

// Add labels to types
// lab def type 1 "Low minimum wage/No change in minw" 	///
//			 2 "Low minimum wage/Change in minw(t+1)"	/// Use this type for leavers specifications
//			 3 "High minimum wage/No change in minw", replace
// lab val type type

// Directories for the public and RDC datasets
*global resultdir 	"results"


// Generate baseline hazard using one of the estimates.
// Note that this

forvalues sex = 1/3 {

// Generate baseline hazard using one of the estimates.
// Note that this

*est use ${resultdir}/joiners_`sex'_1_1_dube_2

// Dube 2 is the specification used to form the baseline hazard.

trial joiners_dube_2 joiners_`sex'_1_1_dube_2


predict xb_bsl,xb // Baseline xb for each bin //

//Next line generates fitted value of minimum wage terms only for each bin.

//Minimum wage terms used in the Dube 2 model.

gen min_bsl =  	   _b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a
				   
// Minimum wage terms with interactions

/*gen min_bsl_int =  	_b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a +
				   	_b[min1_mw]*min1_mw + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1_mw]*min6b1_mw       + _b[min35b1_mw]*min35b1_mw + 	///	
				   _b[min12b1_mw]*min12b1_mw     + _b[min12a1_mw]*min12a1_mw + ///
				   _b[min35a1_mw]*min35a1_mw     + _b[min610a1_mw]*min610a1_mw +  ///
				   _b[min1115a1_mw]*min1115a1_mw + _b[min1620a1_mw]*min1620a1_mw +	///
				   _b[min2125a1_mw]*min2125a1_mw + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1_mw]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1_mw]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin_mw]*cmin_mw +	///
				   _b[cminold2new_mw]*cminold2new_mw  + _b[cminnew_mw]*cminnew_mw + 		///
				   _b[cmin6b_mw]*cmin6b_mw       + _b[cmin35b_mw]*cmin35b_mw + 		///
				   _b[cmin12b_mw]*cmin12b_mw     + _b[cmin12a_mw]*cmin12a_mw + ///
				   _b[cmin35a_mw]*cmin35a_mw     + _b[cmin610a_mw]*cmin610a_mw +  ///
				   _b[cmin1115a_mw]*cmin1115a_mw + _b[cmin1620a_mw]*cmin1620a_mw +	///
				   _b[cmin2125a_mw]*cmin2125a_mw + _b[cmin2630a_mw]*cmin2630a_mw +	///
				   _b[cmin3135a_mw]*cmin3135a_mw + _b[cmin3640a_mw]*cmin3640a_mw + ///
				   _b[cmin4145a_mw]*cmin4145a_mw + _b[cmin4650a_mw]*cmin4650a_mw */
				   
			   

gen xb_latent_bsl = xb_bsl-min_bsl

gen hzd_bsl = 1-exp(-exp(xb_bsl)) //"manual" verification of above.
replace hzd_bsl=. if wagcat==164
gen hzd_latent_bsl = 1-exp(-exp(xb_latent_bsl))

scatter hzd_latent_bsl wagcat if type==1|type==2|type==3
xtset type wagcat

//Smooth baseline hazards

gen     avg_pbin_latent_bsl = 1/3*L1.hzd_latent_bsl + 1/3*hzd_latent_bsl    + 1/3*F1.hzd_latent_bsl
replace avg_pbin_latent_bsl = 3/5*hzd_latent_bsl    + 1/5*F1.hzd_latent_bsl + 1/5*F2.hzd_latent_bsl 					if wagcat ==   1
replace avg_pbin_latent_bsl = 1/5*L1.hzd_latent_bsl + 2/5*hzd_latent_bsl    + 1/5*F1.hzd_latent_bsl + 1/5*F2.hzd_latent_bsl  if wagcat ==   2
replace avg_pbin_latent_bsl = 1/5*F1.hzd_latent_bsl + 2/5*hzd_latent_bsl    + 1/5*L1.hzd_latent_bsl + 1/5*L2.hzd_latent_bsl 	if wagcat == 162
replace avg_pbin_latent_bsl = 3/5*hzd_latent_bsl    + 1/5*L1.hzd_latent_bsl + 1/5*L2.hzd_latent_bsl 					if wagcat == 163
replace avg_pbin_latent_bsl = . if wagcat == 164



replace hzd_latent_bsl=. if wagcat==164

scatter hzd_bsl wagcat if type==1,c(l) msize(vtiny) || scatter hzd_latent_bsl wagcat if type==1,c(l) msize(vtiny) || scatter avg_pbin_latent_bsl wagcat if type==1,c(l) msize(vtiny)


*Can use different models here by adding the suffix in the list.
*Each model type has its own minimum wage terms, which are added to the baseline hazard before forming the fitted densities.

foreach model in spec_4_int_a  {



// Load estimation results for sample 3
*est use ${resultdir}/joiners_`sex'_1_1_`model'
*est store m_`model'

trial Model_estimates_old_interact joiners_`sex'_1_1_`model'

// Minimum wage effects for the current model.

/* These are the original minimum wage-interaction effects - based on levels. */

if "`model'"=="spec_4_int_a" {
gen min_current =  	_b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a + ///
				   _b[min1_mw]*min1_mw + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1_mw]*min6b1_mw       + _b[min35b1_mw]*min35b1_mw + 	///	
				   _b[min12b1_mw]*min12b1_mw     + _b[min12a1_mw]*min12a1_mw + ///
				   _b[min35a1_mw]*min35a1_mw     + _b[min610a1_mw]*min610a1_mw +  ///
				   _b[min1115a1_mw]*min1115a1_mw + _b[min1620a1_mw]*min1620a1_mw +	///
				   _b[min2125a1_mw]*min2125a1_mw + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1_mw]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1_mw]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin_mw]*cmin_mw +	///
				   _b[cminold2new_mw]*cminold2new_mw  + _b[cminnew_mw]*cminnew_mw + 		///
				   _b[cmin6b_mw]*cmin6b_mw       + _b[cmin35b_mw]*cmin35b_mw + 		///
				   _b[cmin12b_mw]*cmin12b_mw     + _b[cmin12a_mw]*cmin12a_mw + ///
				   _b[cmin35a_mw]*cmin35a_mw     + _b[cmin610a_mw]*cmin610a_mw +  ///
				   _b[cmin1115a_mw]*cmin1115a_mw + _b[cmin1620a_mw]*cmin1620a_mw +	///
				   _b[cmin2125a_mw]*cmin2125a_mw + _b[cmin2630a_mw]*cmin2630a_mw +	///
				   _b[cmin3135a_mw]*cmin3135a_mw + _b[cmin3640a_mw]*cmin3640a_mw + ///
				   _b[cmin4145a_mw]*cmin4145a_mw + _b[cmin4650a_mw]*cmin4650a_mw
}

// These are the original minimum wage effects - also commented out.
/*gen min_current =  _b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a */
				   
// Applies current minimum wage effects to the latent baseline model.
		
	
		   
di "here"	

if "`model'"=="v1" {
	gen min_current =  _b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///	   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a
				  
}

if "`model'"=="v1_5" {
	gen min_current =  _b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a 
}


if "`model'"=="1_lag" {		
gen min_current =  _b[min1]*min1 +	///					
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///	   
				   _b[cmin12]*cmin12 +	///
				   _b[cminold2new12]*cminold2new12  + _b[cminold12]*cminold12 + ///
				   _b[cmin6b12]*cmin6b12       + _b[cmin35b12]*cmin35b12 + 		///
				   _b[cmin12b12]*cmin12b12     + _b[cmin12a12]*cmin12a12 + ///
				   _b[cmin35a12]*cmin35a12     + _b[cmin610a12]*cmin610a12 +  ///
				   _b[cmin1115a12]*cmin1115a12 + _b[cmin1620a12]*cmin1620a12 +	///
				   _b[cmin2125a12]*cmin2125a12 + _b[cmin2630a12]*cmin2630a12 +	///
				   _b[cmin3135a12]*cmin3135a12 + _b[cmin3640a12]*cmin3640a12 + ///
				   _b[cmin4145a12]*cmin4145a12 + _b[cmin4650a12]*cmin4650a12 + ///
				   _b[cmin312]*cmin312 +	///
				   _b[cminold2new312]*cminold2new312  + _b[cminold312]*cminold312 + ///
				   _b[cmin6b312]*cmin6b312       + _b[cmin35b312]*cmin35b312 + 		///
				   _b[cmin12b312]*cmin12b312     + _b[cmin12a312]*cmin12a312 + ///
				   _b[cmin35a312]*cmin35a312     + _b[cmin610a312]*cmin610a312 +  ///
				   _b[cmin1115a312]*cmin1115a312 + _b[cmin1620a312]*cmin1620a312 +	///
				   _b[cmin2125a312]*cmin2125a312 + _b[cmin2630a312]*cmin2630a312 +	///
				   _b[cmin3135a312]*cmin3135a312 + _b[cmin3640a312]*cmin3640a312 + ///
				   _b[cmin4145a312]*cmin4145a312 + _b[cmin4650a312]*cmin4650a312
}

if "`model'"=="1_lag_forward" {		
gen min_current =  _b[min1]*min1 +	///					
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///	   
				   _b[cmin12]*cmin12 +	///
				   _b[cminold2new12]*cminold2new12  + _b[cminold12]*cminold12 + ///
				   _b[cmin6b12]*cmin6b12       + _b[cmin35b12]*cmin35b12 + 		///
				   _b[cmin12b12]*cmin12b12     + _b[cmin12a12]*cmin12a12 + ///
				   _b[cmin35a12]*cmin35a12     + _b[cmin610a12]*cmin610a12 +  ///
				   _b[cmin1115a12]*cmin1115a12 + _b[cmin1620a12]*cmin1620a12 +	///
				   _b[cmin2125a12]*cmin2125a12 + _b[cmin2630a12]*cmin2630a12 +	///
				   _b[cmin3135a12]*cmin3135a12 + _b[cmin3640a12]*cmin3640a12 + ///
				   _b[cmin4145a12]*cmin4145a12 + _b[cmin4650a12]*cmin4650a12 + ///
				   _b[cmin312]*cmin312 +	///
				   _b[cminold2new312]*cminold2new312  + _b[cminold312]*cminold312 + ///
				   _b[cmin6b312]*cmin6b312       + _b[cmin35b312]*cmin35b312 + 		///
				   _b[cmin12b312]*cmin12b312     + _b[cmin12a312]*cmin12a312 + ///
				   _b[cmin35a312]*cmin35a312     + _b[cmin610a312]*cmin610a312 +  ///
				   _b[cmin1115a312]*cmin1115a312 + _b[cmin1620a312]*cmin1620a312 +	///
				   _b[cmin2125a312]*cmin2125a312 + _b[cmin2630a312]*cmin2630a312 +	///
				   _b[cmin3135a312]*cmin3135a312 + _b[cmin3640a312]*cmin3640a312 + ///
				   _b[cmin4145a312]*cmin4145a312 + _b[cmin4650a312]*cmin4650a312 + ///
				   _b[cmin6bf12]*cmin6bf12 + _b[cmin35bf12]*cmin35bf12 + ///
				   _b[cmin12bf12]*cmin12bf12 + _b[cminf12]*cminf12 + ///
				   _b[cminold2newf12]*cminold2newf12 + _b[cminnewf12]*cminnewf12 + ///
				   _b[cmin12af12]*cmin12af12 + _b[cmin35af12]*cmin35af12 + ///
				   _b[cmin610af12]*cmin610af12 +_b[cmin1115af12]*cmin1115af12 + ///
				   _b[cmin1620af12]*cmin1620af12 + _b[cmin2125af12]*cmin2125af12 + ///
				   _b[cmin2630af12]*cmin2630af12 + _b[cmin3135af12]*cmin3135af12 + ///
				   _b[cmin3640af12]*cmin3640af12 + _b[cmin4145af12]*cmin4145af12 + ///
				   _b[cmin4650af12]*cmin4650af12	
				   
}

if "`model'"=="seven" {		
gen min_current =  _b[min1]*min1 +	///					
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///	 
				   _b[min5155a1]*min5155a1 + _b[min5660a1]*min5660a1 + ///	
				   _b[min6165a1]*min6165a1 + _b[min6670a1]*min6670a1 + ///	
				   _b[cmin12]*cmin12 +	///
				   _b[cminold2new12]*cminold2new12  + _b[cminold12]*cminold12 + ///
				   _b[cmin6b12]*cmin6b12       + _b[cmin35b12]*cmin35b12 + 		///
				   _b[cmin12b12]*cmin12b12     + _b[cmin12a12]*cmin12a12 + ///
				   _b[cmin35a12]*cmin35a12     + _b[cmin610a12]*cmin610a12 +  ///
				   _b[cmin1115a12]*cmin1115a12 + _b[cmin1620a12]*cmin1620a12 +	///
				   _b[cmin2125a12]*cmin2125a12 + _b[cmin2630a12]*cmin2630a12 +	///
				   _b[cmin3135a12]*cmin3135a12 + _b[cmin3640a12]*cmin3640a12 + ///
				   _b[cmin4145a12]*cmin4145a12 + _b[cmin4650a12]*cmin4650a12 + ///
				   _b[cmin5155a12]*cmin5155a12 + _b[cmin5660a12]*cmin5660a12 + ///
				   _b[cmin6165a12]*cmin6165a12 + _b[cmin6670a12]*cmin6670a12 + ///
				   _b[cmin312]*cmin312 +	///
				   _b[cminold2new312]*cminold2new312  + _b[cminold312]*cminold312 + ///
				   _b[cmin6b312]*cmin6b312       + _b[cmin35b312]*cmin35b312 + 		///
				   _b[cmin12b312]*cmin12b312     + _b[cmin12a312]*cmin12a312 + ///
				   _b[cmin35a312]*cmin35a312     + _b[cmin610a312]*cmin610a312 +  ///
				   _b[cmin1115a312]*cmin1115a312 + _b[cmin1620a312]*cmin1620a312 +	///
				   _b[cmin2125a312]*cmin2125a312 + _b[cmin2630a312]*cmin2630a312 +	///
				   _b[cmin3135a312]*cmin3135a312 + _b[cmin3640a312]*cmin3640a312 + ///
				   _b[cmin4145a312]*cmin4145a312 + _b[cmin4650a312]*cmin4650a312 + ///
				   _b[cmin5155a312]*cmin5155a312 + _b[cmin5660a312]*cmin5660a312 + ///
				   _b[cmin6165a312]*cmin6165a312 + _b[cmin6670a312]*cmin6670a312 + ///
				   _b[cminf12]*cminf12 + ///
				   _b[cminold2newf12]*cminold2newf12 + _b[cminnewf12]*cminnewf12 + ///
				   _b[cmin6bf12]*cmin6bf12 		+ _b[cmin35bf12]*cmin35bf12 + ///
				   _b[cmin12bf12]*cmin12bf12   	+ _b[cmin12af12]*cmin12af12 + ///
				   _b[cmin35af12]*cmin35af12 	+ _b[cmin610af12]*cmin610af12 + ///
				   _b[cmin1115af12]*cmin1115af12 + _b[cmin1620af12]*cmin1620af12 +	///
				   _b[cmin2125af12]*cmin2125af12 + _b[cmin2630af12]*cmin2630af12 +	///
				   _b[cmin3135af12]*cmin3135af12 + _b[cmin3640af12]*cmin3640af12 + ///
				   _b[cmin4145af12]*cmin4145af12 + _b[cmin4650af12]*cmin4650af12 + ///
				   _b[cmin5155af12]*cmin5155af12 + _b[cmin5660af12]*cmin5660af12 + ///
				   _b[cmin6165af12]*cmin6165af12 + _b[cmin6670af12]*cmin6670af12 
			
				   
}

     
if "`model'"=="original" {
gen min_current =  _b[min1]*min1 +	///					
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///	   
				   _b[cmin12]*cmin12 +	///
				   _b[cminold2new12]*cminold2new12  + _b[cminold12]*cminold12 + ///
				   _b[cmin6b12]*cmin6b12       + _b[cmin35b12]*cmin35b12 + 		///
				   _b[cmin12b12]*cmin12b12     + _b[cmin12a12]*cmin12a12 + ///
				   _b[cmin35a12]*cmin35a12     + _b[cmin610a12]*cmin610a12 +  ///
				   _b[cmin1115a12]*cmin1115a12 + _b[cmin1620a12]*cmin1620a12 +	///
				   _b[cmin2125a12]*cmin2125a12 + _b[cmin2630a12]*cmin2630a12 +	///
				   _b[cmin3135a12]*cmin3135a12 + _b[cmin3640a12]*cmin3640a12 + ///
				   _b[cmin4145a12]*cmin4145a12 + _b[cmin4650a12]*cmin4650a12  	
	
}
	

if "`model'"=="extend" {
gen min_current =  _b[min1]*min1 +	///					
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///	
				   _b[min51100a1]*min51100a1 + _b[min101a1]*min101a1 + ///
				   _b[cmin12]*cmin12 +	///
				   _b[cminold2new12]*cminold2new12  + _b[cminold12]*cminold12 + ///
				   _b[cmin6b12]*cmin6b12       + _b[cmin35b12]*cmin35b12 + 		///
				   _b[cmin12b12]*cmin12b12     + _b[cmin12a12]*cmin12a12 + ///
				   _b[cmin35a12]*cmin35a12     + _b[cmin610a12]*cmin610a12 +  ///
				   _b[cmin1115a12]*cmin1115a12 + _b[cmin1620a12]*cmin1620a12 +	///
				   _b[cmin2125a12]*cmin2125a12 + _b[cmin2630a12]*cmin2630a12 +	///
				   _b[cmin3135a12]*cmin3135a12 + _b[cmin3640a12]*cmin3640a12 + ///
				   _b[cmin4145a12]*cmin4145a12 + _b[cmin4650a12]*cmin4650a12  ///	
	
}	
	
gen xb_current=xb_latent_bsl + min_current

/*est use ${resultdir}/joiners_`sex'_1_1_spec_4
est store m_4

gen min_current2 =  _b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a */
				   
*/

scatter min_current wagcat if type==1 & wagcat<=88,c(l) msize(vtiny) || scatter min_current /// 
	wagcat if type==2 & wagcat<=88,c(l) msize(vtiny)|| scatter min_current /// 
	wagcat if type==3 & wagcat<=88,c(l) msize(vtiny)|| scatter min_current /// 
	wagcat if type==4 & wagcat<=88,c(l) msize(vtiny)|| scatter min_current /// 
	wagcat if type==5 & wagcat<=88,c(l) msize(vtiny)


scatter min_current wagcat if type==1 & wagcat<=88,c(l) msize(vtiny) || scatter min_current /// 
	wagcat if type==2 & wagcat<=88,c(l) msize(vtiny) legend(order(1 "$6.75" 2 "7.25 - in 1-2 months" ))
	
	
	
// Probability of being in current bin - this is the hazard using the latent baseline hazard and the 
// Mminimum wage effects of the current model.

gen pbin = 1-exp(-exp(xb_current))
*gen pbin=avg_pbin_latent_bsl * exp(min_current)
 
replace pbin=1 if wagcat==164
gen pbin_approx=exp(xb_current)

table wagcat if type==1,stat(mean pbin) stat(mean pbin_approx)



// Smooth the predicted effect (hazard) of the non minimum wage variables
// within year, province, minwage category, quarter, age group and education group
// (Because we are using types instead of the observed sample, take simply
// the variable type as id variable)
// gen double x_sect = year*1000000000 + prov*10000000 + mincat*100000 + quarter*1000 + agegrp*100 + edgrp*10 + sub

/*gen     avghzdlatu3 = 1/3*L1.hzdlatu3 + 1/3*hzdlatu3    + 1/3*F1.hzdlatu3
replace avghzdlatu3 = 3/5*hzdlatu3    + 1/5*F1.hzdlatu3 + 1/5*F2.hzdlatu3 					if wagcat ==   1
replace avghzdlatu3 = 1/5*L1.hzdlatu3 + 2/5*hzdlatu3    + 1/5*F1.hzdlatu3 + 1/5*F2.hzdlatu3  if wagcat ==   2
replace avghzdlatu3 = 1/5*F1.hzdlatu3 + 2/5*hzdlatu3    + 1/5*L1.hzdlatu3 + 1/5*L2.hzdlatu3 	if wagcat == 162
replace avghzdlatu3 = 3/5*hzdlatu3    + 1/5*L1.hzdlatu3 + 1/5*L2.hzdlatu3 					if wagcat == 163
replace avghzdlatu3 = 1 																	if wagcat == 164
*/
// gen     smthhzdu3   = avghzdlatu3*hzdminu3
// gen smthhzdu3 = avghzdlatu3
// Below are the probabilities (smoothed) of falling in a given bin, conditional 
// on reaching that bin.

// gen avghzdlatu3=hzdlatu3
// gen smthhzdu3=hzdlatu3
// gen pu3 = 1 - exp(-avghzdlatu3)
// gen ps3 = 1 - exp(-smthhzdu3)

// By definition, hazard in top-coded bin is one
// replace pu3=1 if wagcat==164
// replace ps3=1 if wagcat==164
 
// Note pu and pr are hazards
//gen du3   = .
//gen suru3 = .
gen ds3   = .
gen surs3 = .

// For first bin, probs and hzds are the same. 
// replace du3   =     pu3 if wagcat == 1
// replace suru3 = 1 - du3 if wagcat == 1

replace ds3   = pbin if wagcat == 1
replace surs3 = 1 - ds3 if wagcat == 1

// construction of hzds and survivor function for all remaining bins. 

local i = 2
while `i' <= 164 {
	display `i'
	gen probu3    = pbin*L1.surs3
	gen survu3    = L1.surs3 - probu3
	replace ds3   = probu3 if wagcat == `i'
	 
	replace surs3 = survu3 if wagcat == `i'
	drop probu3 survu3	
	 
	local i=`i'+1
}
sum ds3 surs3 if wagcat==164

// Display the results 
// ------------------------------------------------------- //

// There are two ways to present the results. First, we can compare different
// estimation results, conditional on the same types. The models to be compared
// are specified on the global variables $file1 and $file2. Second, compare the
// results for a given sample for different types.

// use $workdatadir\20150115_densities, clear

//gen d  = du1 - ds1
//gen d2 = d*d

// Midpoint is for labeling the x-axis on the wage density polots. 
gen midpoint = .
set more off
replace midpoint = 1.50 if wagcat == 1
replace midpoint = 3.25 if wagcat == 2
replace midpoint = 3.75 if wagcat == 3

local i = 1
while `i' <= 161{
	replace midpoint = 4.05 + (`i'-1)*.1 if wagcat == 3 + `i'
	local i=`i'+1
}

*drop if type==3
table midpoint type,stat(mean ds3)
collect preview
collect export "$dir_path\fitted_densities\joiners_`sex'_1_$j_`model'_`wkdate'.xlsx", open replace sheet(joiners_`sex'_1_$j_`model')


// IMPORTANT: THE NOTES FOR THE GRAPHS HAVE TO BE CHANGED MANUALLY!!!!!

// Compare the estimated CDF 
gen cdfu3 = 1 - surs3	// Based on the survival function, compute the CDF for each type

// gen cdfs3 = 1 - surs3	// Based on the survival function, compute the CDF for each type

// Leavers - densities
/*graph twoway 	(line ds3 midpoint if wagcat <= 163, color(black red) lwidth(0.75 0.75))	///
	if type==1 |type==2 | type==3,																	///
	xtitle("Wage") ytitle("Estimated Density") 											///
	legend(order(1 "No Minimum Wage" 2 "Minimum Wage = 11" ) region(lcolor(none)))		///
	subtitle(, pos(12) nobox) plotregion(lcolor(none))									///
	by(type,																			///
	title("Wage Densities With and Without the Minimum Wage" "Joiners")					///
	note("The selected type of worker is male, 25 years old, with less than High School, living in Ontario"	///
	"during November-2009.")) name(den1, replace) 
	graph export "$dir_path\fitted_densities\joiners_`sex'_1_1_`file'.emf",replace
*/	

line ds3 midpoint if wagcat<=163 & type==1,c(.) || line ds3 midpoint if ///
	wagcat<=163 & type==2,c(.) || line ds3 midpoint if wagcat<=163 & type==3,c(.) ///
	|| line ds3 midpoint if wagcat<=163 & type==4,c(.) ||  line ds3 midpoint if wagcat<=163 & type==5,c(.) ///
	legend(order(1 "$6.75" 2 "7.25 - in 1-2 months" 3 "7.25 - 1-2 months" 4 "7.25 - 3-12 months" 5 "7.25 - >12 months")  ) ///
	xtitle("Wage") ytitle("Estimated Partial Density") 
	graph export "$dir_path\fitted_densities\joiners_`sex'_1_$j_`model'.pdf",replace


gen norm=wagcat-mincat12
replace norm=. if norm<-5 | norm>40 
gen rel=norm*.1
label var rel "Dollars above the minimum wage"
label var min_current "Effect on the probability"
	
graph bar (mean) min_current if type==1, over(norm)  
scatter min_current rel if type==1,c(l) || scatter min_current rel if type==2,c(l) yline(-2(.5)2) ylabel(-2(.5)2) legend(label(1 "Before") label(2 "After"))
*graph export "$dir_path\Fig10a_`model'.eps",replace	
/* CDF plots - commented out
line cdfu3 midpoint if wagcat<=163 & type==1,c(.) || line cdfu3 midpoint if ///
	wagcat<=163 & type==2,c(.) || line cdfu3 midpoint if wagcat<=163 & type==3,c(.) ///
	legend(order(1 "2009 - $11992.78" 2 "1997 - $13821.29" 3 "2002 - $12778.57") ) ///
	xtitle("Wage") ytitle("Estimated Cumulative Density") 
********/	

	drop min_current xb_current pbin pbin_approx midpoint ds3 cdfu3 surs3 norm rel
}



/*_b[min1]*min1 + 									/// Estimated hazard - Effect of minimum wage variables only
				   _b[min6b1]*min6b1       + _b[min35b1]*min35b1 + 	///	
				   _b[min12b1]*min12b1     + _b[min12a1]*min12a1 + ///
				   _b[min35a1]*min35a1     + _b[min610a1]*min610a1 +  ///
				   _b[min1115a1]*min1115a1 + _b[min1620a1]*min1620a1 +	///
				   _b[min2125a1]*min2125a1 + _b[min2630a1]*min2630a1 +	///
				   _b[min3135a1]*min3135a1 + _b[min3640a1]*min3640a1 + ///
				   _b[min4145a1]*min4145a1 + _b[min4650a1]*min4650a1 + ///			   
				   _b[cmin]*cmin +	///
				   _b[cminold2new]*cminold2new  + _b[cminnew]*cminnew + 		///
				   _b[cmin6b]*cmin6b       + _b[cmin35b]*cmin35b + 		///
				   _b[cmin12b]*cmin12b     + _b[cmin12a]*cmin12a + ///
				   _b[cmin35a]*cmin35a     + _b[cmin610a]*cmin610a +  ///
				   _b[cmin1115a]*cmin1115a + _b[cmin1620a]*cmin1620a +	///
				   _b[cmin2125a]*cmin2125a + _b[cmin2630a]*cmin2630a +	///
				   _b[cmin3135a]*cmin3135a + _b[cmin3640a]*cmin3640a + ///
				   _b[cmin4145a]*cmin4145a + _b[cmin4650a]*cmin4650a) */

drop xb_bsl min_bsl hzd_bsl xb_latent_bsl hzd_latent_bsl avg_pbin_latent_bsl

}


log close	
	
/*graph twoway (line hzdminu3  midpoint if wagcat <= 163,color(red) lwidth(0.75 0.75))	 ///
	if type == 1 | type == 2 | type==3 ,									///
	xtitle("Wage") ytitle("Estimated Hazard") 										///
	legend(order(1 "No Minimum Wage" 2 "Minimum Wage=10" ) region(lcolor(none)))	///
	subtitle(, pos(12) nobox) plotregion(lcolor(none))								///
	by(type,																		///
	title("Wage Hazards With and Without the Minimum Wage" "Joiners")				///
	note("The selected type of worker is male, 25 years old, with less than High School, living in Ontario"	///
	"during November-2009.")) name(den3, replace) 
	// graph export ${out_plots}\joiners_baseline_hazard.png,replace	
*/

 /*graph twoway 	(line ds3 midpoint if wagcat <= 163 & type == 1, color(red) lwidth(0.75 0.75))			///
	(line ds3 midpoint if wagcat <= 163 & type == 2, color(black) lwidth(0.75 0.75))			///
	(line ds3 midpoint if wagcat <= 163 & type == 3, color(blue) lwidth(0.75 0.75)),		///
	xtitle("Wage") ytitle("Estimated Density") 														///
	legend(order(1 "Min Wage=10; no changes" 2 "Min Wage change 10 to 11" 3 "Min Wage 11; no change" ) region(lcolor(none)))	///
	title("Wage Densities" "Joiners") plotregion(lcolor(none))		///
	note("The selected type of worker is male, 25 years old, with less than High School, living in Ontario"	///
	"during November-2009.") text(0.03 11 "New W{sub:m} = 11", placement(e)) xline(10.95, lcolor(black) lpattern(dash)) name(dens3, replace)
	graph export ${dir_path}\joiners_baseline_density4lines.png,replace							
*/
/*	
graph twoway 	(line ds3 midpoint if wagcat <= 163 & type == 1, color(red) lwidth(0.75 0.75))			///
	(line ds3 midpoint if wagcat <= 163 & type == 2, color(black) lwidth(0.75 0.75))			///
	(line ds3 midpoint if wagcat <= 163 & type == 3, color(blue) lwidth(0.75 0.75)),		///
	xtitle("Wage") ytitle("Estimated Density") 														///
	legend(order(1 "Min Wage=10; no changes" 2 "Min Wage change 10 to 11" 3 "Min Wage 11; no change" ) region(lcolor(none)))	///
	title("Wage Densities" "Ontario, 2009, ) plotregion(lcolor(none))		///
	note("The selected type of worker is male, 25 years old, with less than High School, living in Ontario"	///
	"during November-2009.")  name(dens3, replace) 
	// graph export ${out_plots}\joiners_baseline_density4lines.png,replace					
		
graph twoway 	(line hzdminu3 midpoint if wagcat <= 163 & type == 1, color(red) lwidth(0.75 0.75))			///
	(line hzdminu3 midpoint if wagcat <= 163 	 & type == 2, color(black) lwidth(0.75 0.75))			///
	(line hzdminu3 midpoint if wagcat <= 163 	 & type == 3, color(blue) lwidth(0.75 0.75)),			///
	xtitle("Wage") ytitle("Estimated Hazard") 														///
	legend(order(1 "Min Wage=10; no changes" 2 "Min Wage change; 10 to 11" 3 "Min Wage 11; no change" ) region(lcolor(none)))	///
	title("Wage Hazards" "Joiners") plotregion(lcolor(none))		///
	note("The selected type of worker is male, 25 years old, with less than High School, living in Ontario"	///
	"during November-2009.") text(6 11 "New W{sub:m} = 11", placement(e)) xline(10.95, lcolor(black) lpattern(dash)) name(hzds3, replace) 
	// graph export ${out_plots}\joiners_baseline_hazard4lines.png,replace	
*/	
	
